% Gerstner's model
clear all

C=281 % *pF
gL=30 %*nS
taum=C/gL
EL=-70.6 %*mV # Same as changing I
VT=-50.4 %*mV
DeltaT=3 %*mV
Vcut=VT+5*DeltaT
b=0 %*nA
Vr=-70.6 %*mV
tauw=taum
a=gL


dt = .1; % ms;

I = 1050
vm(1)=EL;
w(1)=0;
t(1)=0;
for i = 2:ceil(1000/dt)

    dvm=(gL*(EL-vm(i-1))+gL*DeltaT*exp((vm(i-1)-VT)/DeltaT)+I-w(i-1))/C; % volt
    vm(i) = vm(i-1) + dt*dvm;
    dw=(a*(vm(i-1)-EL)-w(i-1))/tauw; % amp
    w(i)=w(i-1) + dt*dw;
    
    if( vm(i) > VT )
        vm(i-1)=0;
        vm(i) = Vr;
        w(i) = w(i)+b;
    end
    t(i)=t(i-1)+dt;
end

plot(t,[vm' w']);
